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Abstract 

We consider the problem of isolating the real roots of a square-free 
polynomial with integer coefficients using (variants of) the continued frac- 
tion algorithm (CP). We introduce a novel way to compute a lower bound 
on the positive real roots of univariate polynomials. This allows us to 
derive a worst case bound of Ob{<1 6 + d 4 r 2 + d 3 r 2 ) for isolating the real 
roots of a polynomial with integer coefficients using the classic variant 
[T] of CF, where d is the degree of the polynomial and r the maximum 
bitsize of its coefficients. This improves the previous bound of Sharma 
[30 j by a factor of d 3 and matches the bound derived by Mehlhorn and 
Ray [19 for another variant of CF; it also matches the worst case bound 
of the subdivision-based solvers. 



1 Introduction 

The problem of isolating the real roots of a square-free polynomial with integer 
coefficients is one of the most well-studied problems in symbolic computation 
and computational mathematics. The goal is to compute intervals with rational 
endpoints that contain one and only one real root of the polynomial, and to have 
one interval for every real root. 

If we restrict ourselves to algorithms that perform computations with ra- 
tional numbers of arbitrary size, then we can distinguish two main categories. 
The first one consists of algorithms that are subdivision-based; their process 
mimics binary search. They bisect an initial interval that contains all the real 
roots until they obtain intervals with one or zero real roots. The different 
variants differ in the way that they count the number of real roots inside an 
interval, for example using Sturm's theorem or Descartes' rule of signs, see also 
Th. [TJ Classical representatives are the algorithms STURM, DESCARTES and 
BERNSTEIN. We refer the reader to [UlQlllSinillllllSlIISI^ and references 
therein for further details. The worst case complexity of all variants in this 
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category is C , s((i 6 + d 4 T 2 ), where d is the degree of the polynomial and r the 
maximum bitsize of its coefficients. Especially, for the STURM solver, recently, 
it was proved that its expected case complexity, if we consider certain random 
polynomials as input, is Ob(t d 2 r), where r is the number of real roots [13]. 
Let us also mention the bitstream version of DESCARTES algorithm, cf. [20] and 
references therein. 

The second category contains algorithms that isolate the real roots of a poly- 
nomial by computing their continued fraction expansion (CF). Since successive 
approximants of a real number define an interval that contains this number, 
CF computes the partial quotients of the roots of the polynomial until the 
corresponding approximants correspond to intervals that isolate the real roots. 
Counting of the real roots is based on Descartes' rule of signs (Th.[l]) and termi- 
nation is guaranteed by Vincent's theorem (Th. [3]). There are several variants 
which they differ in the way that they compute the partial quotients. 

The first formulation of the algorithm is due to Vincent [35] , who computed 
the partial quotients by successive transformations of the form x i-> x + 1. An 
upper bound on the number of partial quotients needed was derived by Us- 
pensky [33] . Unfortunately this approach leads to an exponential complexity 
bound. Akritas pQ, see also [HIE], treated the exponential behavior of CF by 
treating the partial quotients as lower bounds of the positive real roots, and 
computed the bounds using Cauchy's bound. With this approach, c repeated 
operations of the form x h > x + 1 could be replaced by x i— ► x + c. However, 
his analysis assumes an ideal positive lower bound, that is that we can compute 
directly the floor of the smallest positive real root. In [31], it was proven, un- 
der the assumption that Gauss-Kuzmin distribution holds for the real algebraic 
numbers, that the expected complexity of CF is Os(d 4 r 2 ). By spreading the 
roots, the expected complexity becomes Os{d 4 +d 3 r) [32]. The first worst-case 
complexity result of CF, C>B{d 8 T 3 ), is due to Sharma [3D], without any assump- 
tion. He also proposed a variant of CF, that combines continued fractions with 
subdivision, with complexity C>B(d 5 T 2 ). All the variants of CF in [30] compute 
lower bounds on the positive roots using Hong's bound [14], which is assumed 
to have quadratic arithmetic complexity. Mehlhorn and Ray [19] proposed a 
novel way of computing Hong's bound based on incremental convex hull com- 
putations with linear arithmetic complexity. A direct consequence is that they 
reduced the complexity of the variant of CF combined with subdivision [3D] 
to 0(d 4 T 2 ), thus matching the worst case complexity of the subdivision-based 
algorithms. Using [19] and fast Taylor shifts [36], the bound [30] on classical 
variant of CF becomes C*s((i 7 r 3 ). 

As far as the numerical algorithms are concerned, the best known bound for 
the problem is due to Pan [23l [22] and Schonhage [28], see also [29], 6 B {d z T). 
Moreover, it seems that Pan's approach could be improved to C*b(cPt). This 
class of algorithms approximate the roots, real and complex, of the input poly- 
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normals up to a precision. They could be turned to root isolation algorithms by 
requiring them to approximate up to the separation bound, that is the minimum 
distance between the roots. The crux of the algorithms is that they recursively 
split the polynomial until we obtain linear factors that approximate sufficiently 
all the roots, real and complex. We also refer to a recent approach that con- 
centrates only on the real roots [24]. For an implementation of Schonhage's 
algorithm we refer the reader to the routine CPRTS, p. 12 in Addenda, based 
on the multitape Turing machinal. We are not aware of any implementation of 
Pan's algorithm. In the special case where all the roots of the polynomial are 
real, also called the real root problem, dedicated numerical algorithms were pro- 
posed by Reif [25] and Ben-Or and Tiwari [6] for approximating the roots. Nev- 
ertheless, their Boolean complexity is also C*s(d 3 r). Quite recently, Sagraloff 
[27] announced a variant of the bitstream version of DESCARTES algorithm with 
complexity 0B(d 3 T 2 ). 

Our contribution. We present a novel way to compute a lower bound on the 
positive real roots of a univariate polynomial (Lem.[5]). The proposed approach 
computes the floor of the root (possible complex) with the smallest positive real 
part that contributes to the number of the sign variations in the coefficients 
list of the polynomial. Our bound is at least as good as Hong's bound [14] . 
Using this lower bound we improve the worst case bit complexity bound of 
the classical variant of CF, obtained by Sharma [30j, by a factor of d 3 . We 
obtain a bound of 6 B (d 6 + d 4 T 2 ) or O b {N 6 ), where N = max{d,r}, (Th. 0, 
which matches the worst case bound of the subdivision-based solvers and also 
matches the bound due to Mehlhorn and Ray [19| achieved for another variant 
of CF; it also matches the worst case bound of the subdivision-based solvers 

nn na [g ezi eh hsi sni E6] . 

Notation. In what follows Ob, resp. O, means bit, resp. arithmetic, com- 
plexity and the Ob, resp. O, notation means that we are ignoring logarithmic 
factors. For a polynomial A £ deg(A) = d denotes its degree and L (A) = r 
the maximum bitsize of its coefficients, including a bit for the sign. For a E (Q, 
L (a) > 1 is the maximum bitsize of the numerator and the denominator. Let 
M (t) denote the bit complexity of multiplying two integers of size r; using 
FFT, M (r) = 0>b(t). To simplify notation, we will assume throughout the 
paper that lg(deg(^)) = Igd = 0(t) = 0(C(A)). By VAR(A) we denote the 
number of sign variations in the list of coefficients of A. We use A 7 to denote 
the minimum distance between a root 7 of a polynomial A and any other root, 
we call this quantity local separation bound; A = min 7 A 7 is the separation 
bound, that is the minimum distance between all the roots of A. 

J http: //www. iai .uni-boim.de/ ~ schoe/tp/TPpage .html 
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2 A short introduction to continued fractions 



Our presentation follows closely [32j . For additional details we refer the reader 
to, e.g., [37l \7\ [34]. In general, a simple (regular) continued fraction is a 
(possibly infinite) expression of the form 

1 

9o H = [9o,9i. 92, •••], 

9i + — 

92 + • • • 

where the numbers g t are called partial quotients, qt £ Z and 9i > 1 for i > 0. 
Notice that go may have any sign, however, in our real root isolation algorithm 
qo > 0, without loss of generality. By considering the recurrent relations 

P-l = 1, Pq = 90, Pn+l = Qn+l Pn + Pn-1, ^ 
Q-l = 0, Qo = 1, Qn+l = 9n+l Qn + Qn-1, 

it can be shown by induction that R n — = [90, 9i, • • • , 9n], f° r n = 0, 1, 2, . . . . 

If 7 = [9o, gi, • • •] then 7 = 9o + - + ■■■ = q + £~ =1 q^qT and 
since this is a series of decreasing alternating terms it converges to some real 
number 7. A finite section R n = = [90, 9i, • • • , Qn] is called the n— th con- 
vergent (or approximant) of 7 and the tails 7 n +i = [9n+i, 9n+2, • • • ] are known 
as its complete quotients. That is 7 = [go, 9i, • • • , Qn, 7»+i] f° r n = 0,1,2, ... . 
There is an one to one correspondence between the real numbers and the con- 
tinued fractions, where evidently the finite continued fractions correspond to 
rational numbers. 

It is known that Q n > Fn+i and that F n +i < cp n < F n +2, where F n is the 
n— th Fibonacci number and <fi = 1+ 2 V ^ is the golden ratio. Continued fractions 
are the best rational approximation (for a given denominator size). This is as 



follows: g„(Q„ + .+q„) ^ 



1 Qn 



In order to indicate or to emphasize that a partial quotient or an approximant 
belong to a specific real number 7, we use the notation qj and RZ = Pn/Qn, 
respectively. We also use and Rn = Pn /Qn , where k is a non-negative 
integer, to indicate that we refer to the (real part of the) root 7fc of a polynomial 
A. The ordering of the roots is considered with respect to the magnitude of 
their (positive) real part. 



3 Worst case complexity of CF 

Theorem 1 (Descartes' rule of sign). The number R of real roots of A{x) 
in (0, 00) is bounded by VAR(A) and we have R = VAR(A) mod 2. 
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Remark 2. In general Descartes' rule of sign obtains an overestimation of the 
number of the positive real roots. However, if we know that A is hyperbolic, 
i.e. has only real roots, or when the number of sign variations is or 1 then it 
counts exactly. 

The CF algorithm depends on the following theorem, which dates back to 
Vincent's theorem in 1836 [35]. The inverse of Th.[3]can be found in [HI81I21]. 
The version of the theorem that we present is due to Alesina and Galuzzi [5j, 
see also [33l H] HI [2] , and its proof is closely connected to the one and two circle 
theorems (refer to [17, 5j and references therein). 

Theorem 3. [5] Let A^7L\x\ be square-free and let A > be the separation 
bound, i.e. the smallest distance between two (complex) roots of A. Let n be 
the smallest index such that F n -i F n A > where F n is the n-th Fibonacci 
number. Then the map x i-> [cq , C\ , . . . , c n , x] , where Cq , Ci , . . . , c n is an arbitrary 
sequence of positive integers, transforms A{x) toA n {x), whose list of coefficients 
has no more than one sign variation. 

For a polynomial A = Yli-o a i x '> wnere 7 correspond to its (complex) roots, 
the Mahler measure, M(A), of A is M (A) = a d ri[ 7 |>i 111. e -S- EH [37]. 
If we further assume that A G Z[i] and L{A) = r then M (A) < \\A\\ 2 < 
v / d+ l||^||oo = 2Vd+ 1, and so Y\\-y\>i \l\ < 2 T Vd+l. 

We will also use the following aggregate bound. For a proof we refer to 
e.g. [321 El SHI EU [15]. 

Theorem 4. Let A £ Z[jc] such that deg(^4) = d and L (A) = r. Let 7 denotes 
its distinct roots, then 

~\A 7 > 2~ d2 M (Ay 2d o -lgf] A 7 = -^lgA 7 < 3d 2 + 3dlgd+3dr. 

7 77 

3.1 The tree 

The CF algorithm relies on Vincent's theorem (Th. [3]) and Descartes' rule of 
sign (Th. [l]) to isolate the positive real roots of a square-free polynomial A. 
The negative roots are isolated after we perform the transformation x i-» — x; 
hence it suffices to consider only the case of positive real roots throughout the 
analysis . 

The pseudo-code of the classic variant of CF is presented in Alg. [TJ 
Given a polynomial A, we compute the floor of the smallest positive real root 
(plb = Positive Lower Bound). The ideal PLB is a function that can determine 
whether a polynomial has positive real roots, and if there are such roots then 
returns the floor of the smallest positive root of the polynomial. 

Then we perform the transformation x i->- x + b, obtaining a polynomial A^,. 
It holds that VAR(^4) > VAR( J 4 b )- The latter polynomial is transformed to A\ by 
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the transformation x \- > 1 + x and if VhR(Ai) = or VAR( J 4 1 ) = 1, then At, has 
0, resp. 1, real roots greater than 1, or equivalently A has 0, resp. 1, real roots 
greater than b + 1 (Th.[l]). If VAR( J 4 1 ) < VAR( J 4 b ) then (possibly) there are real 
roots of At, in (0, 1), or equivalently, there are real roots of A in (b, b + 1), due 
to Budan's theorem. We apply the transformation x i-> 1/(1 + x) to A^, and we 
get the polynomial A 2 . If VAR^) = or VAR(^4 2 ) = I, A b has 0, resp. 1, real 
root less than 1 (Th. 1), or equivalently A has 0, resp. 1, real root less than 
b + 1, or to be more specific in (b, b + 1) (Th.[l]). If the transformed polynomial, 
Ai and A2, have more than one sign variations, then we apply plb to them and 
we repeat the process. 

Following [U [32j [30] we consider the process of the algorithm as an infinite 
binary tree. The nodes of the tree hold polynomials and (isolating) intervals. 
The root of the tree corresponds to the original polynomial A and the shifted 
polynomial A^,. The branch from a node to a right child corresponds to the 
map x i->- x + 1, which yields polynomial A\, while to the left child to the map 
x t->- 1/(1 + a;), which yields polynomial A 2 . The sequence of transformations 
that we perform is equivalent to the sequence of transformations in Th. 02 and 
so the leaves of the tree hold (transformed) polynomials that have no more than 
one sign variation, if Th. [3] holds. 

A polynomial that corresponds to a leaf of the tree and has one sign vari- 
ation it is produced after a transformation as in Th. \3[ using positive integers 
go, qi, ... , g n . The compact form of this is M : x 1— > g"^^"'^ , where q"~* and 
q 2 - are consecutive convergents of the continued fraction [go , gi , . . . , g n 1. The 
polynomial has one real root in (0,oo), thus the (unordered) endpoints of the 
isolating interval are M(0) = g"~' and M(oo) = 

There are different variants of the algorithm that differ in the way they 
compute PLB. A PLB realization that actually computes exactly the floor of the 
smallest positive real root is called ideal, but unfortunately has a prohibitive 
complexity. 

A crucial observation is that Descartes' rule of sign (Th. [l]), that counts 
the number of sign variations depends not only on positive real roots, but also 
on some complex ones; which have positive real part. Roughly speaking CF 
is trying to isolate the positive real parts of the roots of A that contribute to 
the sign variations. Thus, the ideal PLB suffices to compute the floor of the 
smallest positive real part of the roots of A that contribute to the number of 
sign variations. For this we will use Lem. ]E[ Notice that all the positive real 
roots contribute to the number of sign variation of A, but this is not always the 
case for the complex roots with positive real part. 
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3.2 Computing a partial quotient 



Lemma 5. Let A £ such that deg(A) = d and C(A) = t. We can 

compute the first partial quotient, or in the other words the i?ooj@, c, of the 
real part of the root with the smallest real part, that contributes to the sign 
variations of A in Ob {dr lg c + d 2 lg 2 c) . 

Proof: We compute the corresponding integer using the technique of the expo- 
nential search, see for example [18] . Without loss of generality, we may assume 
that the real root is not in (0, 1), since in this case we should return 0. 

We perform the transformation x m- x + 2° to the polynomial, and then the 
transformation x h-> x + 1. If the number of sign variations of the resulting 
polynomial compared to the original one decreases, then 2° = 1 is the partial 
quotient. If not, then we perform the transformation x i->- x + 2 1 . If the number 
of sign variations does not decrease, then we perform x x + 2 2 . Again if the 
number of sign variations does not decrease, then we perform x i-> x + 2 3 and 
so on. Eventually, for some positive integer k, there would be a loss in the sign 
variations between transformations x H> x + 2 fc_1 and x i-> x + 2 k . In this case 
the partial quotient c, which we want to compute, satisfies 2 fc_1 < c < 2 fc < 2 c. 
The exact value of c is computed by performing binary search in the interval 
[2 fc ,2 fc+1 ]. We deduce that the number of transformations that we need to 
perform is 2k + 0(1) = 21g|_cJ + 0(1). 

In the worst case, each transformation corresponds to an asymptotically fast 
Taylor shift with a number of bitsize (lg c) , which costal B ( M (dr + d 2 lg c) lg d) 
[36, Th. 2.4]. By considering fast multiplication algorithms the costs becomes 
Ob {dr + d 2 lg c) and multiplying by the number of transformations needed, lg c, 
we conclude the proof. 

It is worth noticing that we do not consider the cases c = 2 fc or c = 2 fc+1 , 
since then we have computed, exactly, a rational root. □ 

3.3 Shifts operations and total complexity 

Up to some constant factors, we can replace A in Th.[3]by A 7 , refer to [30] for a 
proof. This allows us to estimate the number, m 7 , of partial quotients needed, 
in the worst case, to isolate the positive real part of a root 7. It holds 

m 7 < i(l + log^2-lgA 7 ) < 2 - ^lgA 7 . 

2 We choose to use c instead of go because in the complexity analysis that follow A could be a 
result of a shift operation, thus the computed integer may not be the 0-th partial quotient of the 
root that we are trying to approximate. 

3 Following Th. 2.4(E) in | 36| the cost of performing the operation f(x + a), where deg(/) = n, 
£ (/) = r and L (a) = cr is Os(M {jit + n 2 <r) lg ), and if we assume fast multiplication algorithms 
between integers, then it becomes Og(nT + n 2 cr). 
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The transformed polynomial has either one or zero sign variation and if 7 G 
M, then the corresponding interval isolates 7 from the other roots of A. The 
associated continued fraction of (the real part of) 7 is [q^ , qj, • ■ • , <Jm 7 ]- It holds 
that J2 7 m f = 0{d 2 + dr) [32l[30]. The following lemma bounds the bitsize of 
the partial quotients, qj,, of a root 7. 

Lemma 6. Let A £ such that deg(A) = d and C(A) = r. For the real 

part of any root 7 it holds 

E^J) = ig(^) + E^i) ^ igc^) + 1 - ig a 7 , 



where we assume that > 0, and the term 1 — lg A 7 appears only when A 7 < 1, 
over a subset of distinct roots of A, then 



i.e. when m 1 > 1. Moreover lg(<?o) ^ lg 11^112 < t + lgd and if 7 ranges 



EE's^ < i+r+igd-i g nA 7 =o(<i 2 +dT). 

7 fc=0 7 

Proof: The Mahler measure, M{A), of A is (.4) = a d ri[ 7 |>i M- Jt also 
holds < ||^[1 2 < % /dTT|| J 4|| 00 = 2 r V'dTT, and so It| < 2 T % /dTT. 

Since g 7 is the integer part of 7 it holds Y[ 7 3o < l~I| 7 |>i I T I < ll^lb and thus 

E^o) <lg%/^TT + Ig||i4||oo <r + lgd. (2) 

7 

Following [30j we know that 

^ryyi — > «■ QL,Ql y ~i < 2/a 7 . (3) 

From Eq. we get Q fc = qkQk-i + Qk-2 => Qfc > qkQk-i, for ft > 1. I we 
apply the previous relation recursively we get n*£i 3fc ^ Qm 7 < 2/A 7 and 
llfci 1 ll < <&,-! < 2/A 7 , and so 



Eig^igEK^ 1 - 1 ^- 

fc=l fc=l 
Finally, we sum over all roots 7 and we use ([2]) and Th. HI 

E E is «z = E is so 7 + E E is «z < E * + E ( 1 - * A . ) 

7 fc=0 7 7 fc=l 7 7 

< 1 + t + lg d + d 2 + 3dlg d + 3dr, 
which completes the proof. □ 
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At each step of CF we compute a partial quotient and we apply a Taylor shift 
to the polynomial with this number. In the worst case we increase the bitsize 
of the polynomial by an additive factor of 0(dlg(q1)), at each step. The overall 
complexity of CF is dominated by the computation of the partial quotients. 

The following table summarizes the costs of computing the partial quotients 
of 7 that we need: 

th step 6 B {dT\g{ql) + d 2 \g(q2) lg(<£)) 

1 st step 6 B (dr lg(g7) + d 2 lg(g V) lg(??)) 

(= 6 B (d(r + dlg(tf )) lg(qj) + d 2 lg 2 (g7))) 

2 nd step 6 B {dT \g{ql) + d 2 \g{q2qjql) \g{ql)) 



step d B {dr\g{ql l ) + d 2 \g{Y\ql)\g{q'L)) 



fc=0 



We sum over all steps to derive the cost for isolating 7, C 7 , and after applying 
some obvious simplifications and use Lem.[6]we get 

dr Y, Ig(gZ) + d 2 E Ig EI «7 = Q b [dr £ ^) + W £ ^fa 

fc=0 fc=0 j=0 J \ k=0 \fc=0 

= 5 B (dr(lg(q2) - lg A 7 ) + d 2 (lg 2 (g 7 ) + lg 2 A 7 )) . 

To derive the overall complexity, C, we sum over all the roots that CF tries 
to isolate and we use Lem. M and Th. ffl Then 



c = y CT 

= O s (dr lg(<j 7 ) - dr ]T 7 lg A 7 + d 2 £ 7 lg 2 foj) + d 2 £ 7 lg 2 A 7 ^ 
= O s (drE 7 lg( g J)-drJ] 7 lgA 7 + d 2 (j: 7 lg( g 7 )) 2 +d 2 (j: 7 lgA 7 ) 2 ) ( 4 ) 
= Ob (dr(r + lgd) + dr(d 2 + dlgd + dr) + d 2 (r 2 + lg 2 d) + d 2 (d 4 + d 2 T 2 )) 

= B (d 6 + d 4 T 2 ). 



In the previous equation it possible to write Yl 7 lg ^7 ^ (S 7 lg ^ 7 
because A 7 < 1, and hence lgA 7 < 0, for all 7 that are involved in the sum. 
For the roots that holds A 7 > 1 the algorithm isolates them without computing 
any of their partial quotients, with the exception of g 7 . 

The previous discussion leads to the following theorem. 

Theorem 7. Let A £ 7L [x], where deg(j4) = d and C {A) = r. The worst case 
complexity of isolating the real roots of A using the CF is <D B (d 6 + d 4 T 2 ). 
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Algorithm 1: CF(A, M) 



Input: A G 7L[X],M{X) = fc.I.m.n G Z 

Output: A list of isolating intervals 

Data: Initially M(Jf) = X, i.e. fc = n = 1 and I = m = 

1 if A(0) = then 

2 OUTPUT Interval( M(0), M(0)) ; 

3 A<-A(X)/X; 

4 CF(A,M); 

5 V <- Var(^); 

6 if V = then RETURN; 
r if V = 1 then 

8 OUTPUT Interval( M(0), M(oo)); 

9 RETURN; 

10 b <- PLB(^4) // PLB = PositiveLowerBound ; 
n if b > 1 then A b <- A{b + X), M <- M(b + X) ; 

12 Ax <r- A b {\ + X), Mi <- M(l + X) ; 

13 CF(j4!,Mi) // Looking for real roots in (l,+oo); 

14 A 2 «- A(i^), M 2 <- M( T ^ T ) ; 

15 CF( J 4 2 ,M 2 ) // Looking for real roots in (0,1) ; 

16 return; 
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